#include <src/pcps.h>
#include <src/input.h>
#include <src/hamiltonian.h>
#include <src/wavefunction.h>
#include <src/correlator.h>
#include <src/cluster.h>
#include <src/expectation_value.h>
#include <src/parallel.h>

template <class SCALAR>
SCALAR pcps::energy(const pcps::Input & userinp,
                    const pcps::Hamiltonian<SCALAR> & ham,
                    const pcps::Wavefunction<SCALAR> & wfn,
                    const int perturbation_level)
{

  // define the worker class that will compute the energy in parallel
  class Local_Worker : public pcps::Parallel_Worker_Base<SCALAR> {

    public:
      void do_work(const pcps::Input & userinp,
                   const pcps::Hamiltonian<SCALAR> & ham,
                   const pcps::Wavefunction<SCALAR> & wfn,
                   const int perturbation_level,
                   const boost::shared_array<SCALAR> & output_array,
                   const int array_size,
                   const int my_min,
                   const int my_max,
                   const std::vector<int> & slaves,
                   int & slaves_used,
                   int & jobs_in_progress)
      {
        // compute each operator's contribution to the energy
        for (int i = my_min; i < my_max; i++)
          pcps::expectation_value<SCALAR>(slaves, slaves_used, jobs_in_progress, output_array.get(), 0, pcps::unity<SCALAR>(),
                                          userinp, wfn, ham.ops().at(i), perturbation_level,
                                          pcps::Cluster(), pcps::Cluster(), typename pcps::Correlator<SCALAR>::Set());
      }

  };

  // create the worker
  Local_Worker worker;

  // compute the energy in parallel
  boost::shared_array<SCALAR> energy_array = pcps::compute_in_parallel(userinp, ham, wfn, perturbation_level, &worker, 1);

  // get MPI info
  const MPI::Comm & comm = MPI::COMM_WORLD;
  const int nproc = comm.Get_size();
  const int myrank = comm.Get_rank();

  // broadcast the energy to all processes
  SCALAR retval = ( myrank == 0 ? energy_array[0] : pcps::zero<SCALAR>() );
  comm.Barrier();
  comm.Bcast((void *)&retval, sizeof(SCALAR), MPI::CHAR, 0);

  // return the energy
  return retval;

}

// instantiate the template
template double               pcps::energy(const pcps::Input & userinp,
                                           const pcps::Hamiltonian<double> & ham,
                                           const pcps::Wavefunction<double> & wfn,
                                           const int perturbation_level);
template std::complex<double> pcps::energy(const pcps::Input & userinp,
                                           const pcps::Hamiltonian<std::complex<double> > & ham,
                                           const pcps::Wavefunction<std::complex<double> > & wfn,
                                           const int perturbation_level);
